
program agent_based_termite

  use ifport

  implicit none
  logical :: show_termites=.true., slow=.true.
  integer, parameter :: tslow=1
  real*8, parameter :: dim_c=1d0, sst=dim_c/5
  real*8, parameter :: L=15d0/2 !External arena diameter 
  real*8, parameter :: L2=L*L
  real*8, parameter :: Lin=6.25d0/2 !Internal arena diameter
  integer, parameter :: Nta=1 !Number of runs (samples)
  integer, parameter :: Nc=12 !Number of termites
  integer, parameter :: Nob=1
  integer, parameter :: kp1=1, kp2=Nc
  real*8, parameter :: rho=1.d0*Nc/(4*L2)
  real*8, parameter :: dtp=1e0, dtp_win=dtp
  integer, parameter :: tmax=2e6, Nump_max=1e4, Np=Nc !int(4*rho*L2)
  real*8, parameter :: alpha1=-2.60d0, gamma1=1d0/(1+alpha1), wtmax=tmax
  real*8, parameter :: sigma=0.1d0, ptheta=1e-6, pturn=0.1d0
  real*8, parameter :: alpha=1.8d0, gamma=1d0/(1+alpha)
  real*8, parameter :: pcrow=0.d0, pwt=0d0, step_t=0.0026d0, dt=1d0/Np
  real*8, parameter :: sigma_step=0.2d0, step_mean=0.15d0, p_recuo=0.d0

  integer*1 :: gera, gera1
  integer :: i, win1, win2, win3, k, k1, Nts, Na, Nump, Ncsav
  integer :: iseed, key
  integer*4 :: delay
  real :: xs, ys, rt
  real*8 :: wt(Np), ri(Np,2), ri_old(2,2)!, r1(Nump_max,Nob), r2(Nump_max,Nob), ts(Nump_max)
  real*8 :: dx(Nob), dy(Nob), x, y, r(Nob)
  real*8 :: theta(Np), pi, z, z1, phi, thetan, zgauss, doispi, zgauss_salto
  real*8 :: t, tp, tp1, cx, sy, dx1, dy1, sx, cy
  real*8 :: c1, c2, sinal, ttal, salto, dist, step
  logical :: sobrepos
  character :: str1*8

  rt=dim_c/2
  Ncsav=min(10,Nc)
  delay=0
  iseed=-14571
  z=ran2(iseed)
  do i=1,0
     z=ran2(iseed)
  enddo
  pi=atan(1d0)*4
  doispi=pi*2
  c1=wtmax**(1+alpha1)-1

  call gopen(800,800,win1)
  call gopen(350,350,win2)
  call gopen(350,350,win3)
  xs=-L-2; ys=L+2
  call newwindow(win1,xs,xs,ys,ys)
  call newwindow(win2,xs,xs,ys,ys)
  call newwindow(win3,xs,xs,ys,ys)
  call gsetnonblock(1)
  do Na=1,Nta
     print*, 'Amostra',Na
     call gclr(win1); call gclr(win2); call gclr(win3)
     xs=L+0.2
     call newpencolor(win2,3)
     call drawcirc(win2,0.0,0.0,xs,xs)
     call newpencolor(win3,3)
     call drawcirc(win3,0.0,0.0,xs,xs)
     call init
     call evolui
  enddo
contains
  !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  subroutine init
    real*8 :: dx, dy, dx2, dy2, d, dx1, dy1
    ri=0
    do k=1,Np
1      z=ran2(iseed); z=z*(L-Lin)+Lin
       phi=(ran2(iseed)-0.5)*doispi
       x=z*cos(phi); y=z*sin(phi)
       dx=x*x; dy=y*y
       d=sqrt(dx+dy)
       if (d>L) goto 1
       sobrepos=.false.
       do i=1,k-1
          dx=ri(i,1)-x; dx2=dx*dx
          dy=ri(i,2)-y; dy2=dy*dy
          d=dx2+dy2
          if (d<dim_c) then
             sobrepos=.true.
             exit
          endif
       enddo
       if (sobrepos) goto 1
       ri(k,1)=x
       ri(k,2)=y
       z=ran2(iseed)-0.5; z=z*2
       theta(k)=z*doispi
       z=ran2(iseed)
       wt(k)=-1
       !From active to inactive state
       z1=ran2(iseed)
       if (z1<pwt) then
          z1=ran2(iseed); z1=(c1*z1+1)**gamma1
          wt(k)=z1
       endif
    enddo
    ri_old(1,1)=ri(kp1,1); ri_old(1,2)=ri(kp1,2)
    ri_old(2,1)=ri(kp2,1); ri_old(2,2)=ri(kp2,2)
    call update_win1
  end subroutine init
  !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  subroutine evolui
    integer :: img, tt, nstat
    real*8 :: r, dx, dy, phiL, tp2
    character :: str*5
    t=0; tp=dtp; tp1=0; tt=100; write(str1,'(1I8)')tt
    tp2=0
    gera=1; gera1=1; Nts=0; Nump=0
    img=0
    k=-1
    do_evolui:do
       if (t>tp) then
          tp=tp+dtp;
          if (t>tp1) then
             call ggetch(key)
             if (key/=-1) exit
             tp1=tp1+dtp_win
             call gclr(win1)
             str1=trim(adjustl(str1))
             call drawstr(win1,-20,-20,16.0,str1,0.1,len(str1))
             if (show_termites) call update_win1
          endif
       endif
       do k1=1,Np
          if (wt(k1)>0) wt(k1)=wt(k1)-dt
       enddo
       if (t>tt) then
          write(str1,'(1I8)')tt
          tt=tt+100
       endif
       t=t+dt
       if (t>tmax) exit
       z1=ran2(iseed)
       k=min(Nc,int(z1*Nc+1)) !Random selection of an individual
       !Waiting time
       if (wt(k)>0) then
          !Inactive state
          cycle do_evolui
       else
          !From active to inactive state
          z1=ran2(iseed)
          if (z1<pwt) then
             z1=ran2(iseed); z1=(c1*z1+1)**gamma1
             wt(k)=z1
             cycle do_evolui
          endif
       endif

       !Turning angle determination
       call turn_angle
       theta(k)=theta(k)+z
       !Walk step
       salto=gaussian_fly()

       step=min(sst,salto)
       cx=step*cos(theta(k)); sy=step*sin(theta(k))
       x=ri(k,1); y=ri(k,2); dist=0
       exec_salto:do
2         x=x+cx; y=y+sy
          xs=x; ys=y
          dist=dist+step
          if (dist>=salto) then
             dist=dist-step
             x=x-cx; y=y-sy
             exit exec_salto
          endif
          r=sqrt(x*x+y*y)
          if (r>L) then
             x=x-cx; y=y-sy
             sinal=1
             if (y>=0) then
                thetan=atan2(y,x)
             else
                thetan=doispi+atan2(y,x)
             endif
             sinal=-sign(sinal,cx*sin(thetan)-sy*cos(thetan))
             call turn_angle
             theta(k)=theta(k)+sinal*abs(z)
             cx=step*cos(theta(k)); sy=step*sin(theta(k))
             cycle exec_salto
          else if (r<Lin) then
             x=x-cx; y=y-sy
             sinal=1
             if (y>=0) then
                thetan=atan2(y,x)
             else
                thetan=doispi+atan2(y,x)
             endif
             sinal=sign(sinal,cx*sin(thetan)-sy*cos(thetan))
             call turn_angle
             theta(k)=theta(k)+sinal*abs(z)
             cx=step*cos(theta(k)); sy=step*sin(theta(k))
             cycle exec_salto
          endif
          call verify_sobrepos
          if (sobrepos) then
             cycle do_evolui
             x=x-cx; y=y-sy
             z1=ran2(iseed)
             if (z1<pturn) then
                call turn_angle
                theta(k)=theta(k1)+z
                cx=step*cos(theta(k)); sy=step*sin(theta(k))
                cycle !exit exec_salto
             endif
             z1=ran2(iseed)
             if (z1<pcrow) then
                z1=ran2(iseed); z1=(c1*z1+1)**gamma1
                wt(k1)=z1
                z1=ran2(iseed); z1=(c1*z1+1)**gamma1
                wt(k)=z1
             endif
             exit exec_salto
          endif
          t=t+dt
          !Uncomment this block to write the termite location, termite 1 in fort.11, termite 2 in fort.12 and so on
          !if (t>tp2) then
          !   tp2=tp2+0.1d0
          !   do k1=1,Nc
          !      write(10+k1+(na-1)*Nc,*)t,ri(k1,1),ri(k1,2)
          !   enddo
          !endif
      enddo exec_salto
       ri(k,1)=x; ri(k,2)=y
       if (k==kp1) then
          call newpencolor(win2,2)
          xs=ri_old(1,1); ys=ri_old(1,2)
          call moveto(win2,xs,ys)
          xs=ri(k,1); ys=ri(k,2)
          call lineto(win2,xs,ys)
       else if (k==kp2) then
          call newpencolor(win3,4)
          xs=ri_old(2,1); ys=ri_old(2,2)
          call moveto(win3,xs,ys)
          xs=ri(k,1); ys=ri(k,2)
          call lineto(win3,xs,ys)
       endif
       if (k==kp1) then
          ri_old(1,1)=ri(k,1); ri_old(1,2)=ri(k,2)
       else if (k==kp2) then
          ri_old(2,1)=ri(k,1); ri_old(2,2)=ri(k,2)
       endif
    enddo do_evolui

  end subroutine evolui
  !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  subroutine turn_angle
    real*8 :: r, phi
    z1=ran2(iseed)
    if (z1<ptheta) then
       !Random Pertubation
       z=ran2(iseed)-0.5d0; z=z*doispi
    else
       !Gaussian distribution
       if (gera/=0) then
          z=ran2(iseed); z1=ran2(iseed)
          r=sqrt(-2.*log(z))
          phi=pi*2*z1
          z=r*sin(phi)*sigma
          zgauss=r*cos(phi)*sigma
          gera=0
       else
          z=zgauss
          gera=1
       endif
    endif
  end subroutine turn_angle
  !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  function gaussian_fly()
    real*8 :: gaussian_fly
    real*8 :: theta, r
    if (.not.gera1) then
       gera1=.true.
       gaussian_fly=zgauss_salto
       if (gaussian_fly<0) goto 1
       return
    endif
1   z1=ran2(iseed)
    zgauss_salto=ran2(iseed)
    r=sqrt(-2.*log(z1))
    theta=doispi*zgauss_salto
    z1=r*sin(theta)*sigma_step+step_mean
    zgauss_salto=r*cos(theta)*sigma_step+step_mean
    gaussian_fly=z1
    if (zgauss_salto<0) goto 1
    gera1=.false.
  end function gaussian_fly
  !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  subroutine update_win1
    integer :: k
    call newpencolor(win1,10)
    xs=L+0.55
    call drawcirc(win1,0.0,0.0,xs,xs)
    xs=Lin-0.55
    call drawcirc(win1,0.0,0.0,xs,xs)
    call newlinewidth(win1,1)
    do k=1,Np
       xs=ri(k,1); ys=ri(k,2)
       call newpencolor(win1,1)
       if (wt(k)<0) then
          call newpencolor(win1,11)
       else
          call newpencolor(win1,9)
       endif
       if (k==kp1) then
          call newpencolor(win1,2)
       else if (k==kp2) then
          call newpencolor(win1,4)
       endif
       call draw_termite(k)
    enddo
    if (slow) call msleep(tslow)
    call newlinewidth(win1,1)
  end subroutine update_win1
  !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  subroutine draw_termite(k)
    integer :: k
    real :: dx1, dy1, zv1, zv2
    character :: str*3
    
    if (wt(k)<0) then
       zv1=(ran2(iseed)*1.0+0.0); zv2=sqrt(1-zv1*zv1)
    else
       zv1=0.5; zv2=sqrt(1.0-zv1*zv1)
    endif
    dx1=0.25d0*cos(theta(k)); dy1=0.25d0*sin(theta(k))
    call drawarrow(win1,xs,ys,xs+dx1,ys+dy1,-zv1,zv2,114)
    dx1=0.3d0*cos(theta(k)); dy1=0.3d0*sin(theta(k))
    call drawarrow(win1,xs-dx1/1.5,ys-dy1/1.5,xs+dx1*0,ys+dy1*0,-zv1*2,zv2*2,114)
    if (wt(k)<0) then
       zv1=(ran2(iseed)*0.2+0); zv2=sqrt(0.6-zv1*zv1)
    else
       zv1=0.1; zv2=sqrt(0.6-zv1*zv1)
    endif
    call drawarrow(win1,xs-dx1/1.5,ys-dy1/1.5,xs-dx1*0,ys-dy1*0,-zv1*3,zv2*3,114)
    if (wt(k)<0) then
       zv1=(ran2(iseed)*0.4+0.1); zv2=sqrt(0.6-zv1*zv1)
    else
       zv1=0.3; zv2=sqrt(0.6-zv1*zv1)
    endif
    call drawarrow(win1,xs-dx1/1.5,ys-dy1/1.5,xs-dx1*0,ys-dy1*0,zv1*3,zv2*3,114)
    call fillcirc(win1,xs,ys,0.10,0.10)
    dx1=0.1d0*cos(theta(k)); dy1=0.1d0*sin(theta(k))
    call fillcirc(win1,xs-dx1,ys-dy1,0.12,0.12)
    dx1=0.2d0*cos(theta(k)); dy1=0.2d0*sin(theta(k))
    call fillcirc(win1,xs-dx1,ys-dy1,0.12,0.12)
    dx1=0.2d0*cos(theta(k)); dy1=0.2d0*sin(theta(k))
    call fillcirc(win1,xs+dx1,ys+dy1,0.08,0.08)
    dx1=0.10d0*cos(theta(k)); dy1=0.10d0*sin(theta(k))
    call fillcirc(win1,xs+dx1,ys+dy1,0.05,0.05)
    !return
    xs=ri(k,1); ys=ri(k,2)
    call newpencolor(win1,1)
    call drawcirc(win1,xs,ys,rt,rt)
    !zv1=0.5; zv2=sqrt(1.0-zv1*zv1)
    !dx1=0.45d0*cos(theta(k)); dy1=0.45d0*sin(theta(k))
    !call drawarrow(win1,xs-dx1,ys-dy1,xs+dx1,ys+dy1,0.3,0.4,114)
    write(str,'(1I3)')k
    call drawstr(win1,xs,ys,16.0,str,0.0,4)
  end subroutine draw_termite
  !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  subroutine verify_sobrepos
    integer :: i
    real*8 :: dx, dy, dx2, dy2, d
    sobrepos=.false.
    do i=1,Np
       if (i/=k) then
          dx=ri(i,1)-x; dx2=dx*dx
          dy=ri(i,2)-y; dy2=dy*dy
          d=sqrt(dx2+dy2)
          if (d<dim_c) then
             sobrepos=.true.
             k1=i
             return
          endif
       endif
    enddo
  end subroutine verify_sobrepos
  !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  subroutine verify_sobrepos1
    integer :: i
    real*8 :: dx, dy, dx2, dy2, d
    sobrepos=.false.
    do i=1,Np
       if (i/=k1) then
          dx=ri(i,1)-ri(k1,1); dx2=dx*dx
          dy=ri(i,2)-ri(k1,2); dy2=dy*dy
          d=dx2+dy2
          if (d<dim_c) then
             sobrepos=.true.
             return
          endif
       endif
    enddo
  end subroutine verify_sobrepos1
  !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  FUNCTION ran2(idum)
    INTEGER :: idum
    REAL*8 ::ran2
    Integer, PARAMETER ::IM1=2147483563,IM2=2147483399,IMM1=IM1-1,&
         IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,&
         NTAB=32,NDIV=1+IMM1/NTAB
    real*8,parameter   ::EPS=1.2e-7,RNMX=1.-EPS,AM=1./IM1
    INTEGER            ::idum2,j,k,iv(NTAB),iy
    SAVE iv,iy,idum2
    DATA idum2/123456789/, iv/NTAB*0/, iy/0/
    if (idum.le.0) then
       idum=max(-idum,1)
       idum2=idum
       do j=NTAB+8,1,-1
          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1
          if (idum.lt.0) idum=idum+IM1
          if (j.le.NTAB) iv(j)=idum
       enddo
       iy=iv(1)
    endif
    k=idum/IQ1
    idum=IA1*(idum-k*IQ1)-k*IR1
    if (idum.lt.0) idum=idum+IM1
    k=idum2/IQ2
    idum2=IA2*(idum2-k*IQ2)-k*IR2
    if (idum2.lt.0) idum2=idum2+IM2
    j=1+iy/NDIV
    iy=iv(j)-idum2
    iv(j)=idum
    if(iy.lt.1)iy=iy+IMM1
    ran2=min(AM*iy,RNMX)
    return
  END Function ran2
end program agent_based_termite
